Forced Stratified Turbulence: 
Successive Transitions with Reynolds Number 
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Numerical simulations are made for forced turbulence at a sequence of increasing values of 
Reynolds number, R, keeping fixed a strongly stable, volume-mean density stratification. At smaller 
values of R, the turbulent velocity is mainly horizontal, and the momentum balance is approximately 
cyclostrophic and hydrostatic. This is a regime dominated by so-called pancake vortices, with only 
a weak excitation of internal gravity waves and large values of the local Richardson number, Ri, 
everywhere. At higher values of R there are successive transitions to (a) overturning motions with 
local reversals in the density stratification and small or negative values of Ri; (b) growth of a 
horizontally uniform vertical shear flow component; and (c) growth of a large-scale vertical flow 
component. Throughout these transitions, pancake vortices continue to dominate the large-scale 
part of the turbulence, and the gravity wave component remains weak except at small scales. 

PACS numbers: 47.20.-k, 47.27.-i, 47.55. Hd 
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I. INTRODUCTION 



Most atmospheric and oceanic flows on intermediate 
scales of 0(1) — O(10 3 ) m are strongly influenced by the 
stable vertical density stratification, dp(z)/dz < (with 
z = e 3 anti-parallel to gravity), but are less influenced 
by Earth's rotation than flows on larger scales. Con- 
sequently, the turbulence in this regime is quite differ- 
ent from unstratified shear turbulence, convection, and 
geostrophic turbulence, all of which have been more ex- 
tensively studied and whose behaviors are now relatively 
familiar. It has been argued 13] that the main effect of 
strongly stable stratification — i.e., small Froude num- 
ber, Fr — V/NL V , where V is a horizontal velocity scale, 
N = —(g I ' p ) dp/dz is the Brunt Vaissala frequency, and 
L v is vertical length scale — is to organize the flow into 
two distinct, non-interacting classes: nearly linear inter- 
nal gravity waves and fully nonlinear stratified turbu- 
lence. The flow patterns of stratified turbulence are of- 
ten called pancake vortices 01 or vortical modes 
because of their small aspect ratio, H/L (where L is a 
horizontal length scale), and significant vertical vortic- 
ity, neither of which is generally true for internal grav- 
ity waves. Pancake vortices have an anisotropic veloc- 
ity (primarily horizontal) and shear field (primarily ver- 
tical), and they evolve principally under the nonlinear 
influence of horizontal advection as in two-dimensional 
turbulence. These motions cause little vertical turbulent 
heat or mass flux, and they have a highly anisotropic, 
inhomogeneous energy cascade to small scales and dissi- 
pation 7J. At moderate values of Reynolds number - 
R = VLh/v, where v is the kinematic viscosity and Lh is 
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a horizontal length scale — stratified turbulence evolves 
self-consistently, at least in freely decaying flows, in the 
sense that a bulk value for Fr does not increase with 
time as energy dissipation causes R to decrease Eft At 
leading order in Fr, the inviscid dynamical balances for 
stratified turbulence are equivalent to two-dimensional 
turbulence |l3j evolving independently at each vertical 
level. The energy dissipation may be modeled by adding 
a vertical eddy diffusion Q that acts to couple vertically 
adjacent levels and diffusively selects a limiting vertical 
length scale. However, with the assumption of uniform 
asymptotic validity as Fr — > 0, stratified turbulence is 
constrained by hydrostatic and cyclostrophic force bal- 
ances that also act to couple adjacent layers and may 
internally select a finite vertical scale as R — > 00 without 
inducing any vertically overturning motions [uH ]. 

These behaviors have been studied both with labora- 
tory experiments 0, l24j and with numerical simulations 
up to R = O(100); a review of this subject has recently 
been made by Riley & Lelong ^3- In the ocean and at- 
mosphere, R values are generally several orders of magni- 
tude larger than those commonly reached in experiments 
or numerical simulations. Thus, an important open ques- 
tion is whether the preceding wave-turbulence partition 
remains valid at very large R. A central part of this ques- 
tion is whether the pancake vortices persist and remain 
"stable" with respect to overturning motions. 

The dynamical stability properties of a stably stratified 
shear flow are usually related to the Richardson number, 
Ri, locally defined by 



Ri = 



Po dz 



dz 



= N + — 



duh 

dz 



, (1) 

where u/j is the horizontal velocity and 9 = —gp'/p is 
the normalized "temperature" associated with density 
fluctuations p'. Alternatively, Ri is defined as a bulk 
measure in an analogous fashion with N 2 in the numer- 
ator and shear variance in the denominator (i.e., a bulk 
Ri ~ Fr- 2 ). In the inviscid limit, a sufficient condition 
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for stability of a parallel vertically sheared flow (i.e., a 
Kelvin-Helmholtz flow) is that the local Ri exceed 0.25 
throughout the flow 0,01 . Gage Q refined this criterion 
for several simple viscid shear flows and obtained values 
of the critical Ri for linear instability between about 0.05 
and 0.11 for large R. In a numerical simulation at very 
high resolution (i.e., with a maximum Reynolds number 
based on the shear layer thickness Ren = 24000) and 
moderate stratification, Werne & Fritts [23| show that 
the turbulence organizes itself so that Ri never exceeds 
0.25. In the more complex natural environment, velocity 
variance and estimates of the vertical mixing efficiency 
increase rapidly as Ri decreases through the range be- 
tween about 0.5 and 0.0 0. On the other hand, Majda 
& Shefter stress the importance of temporal behavior 
on flow stability by constructing a family of time periodic 
solutions that are unstable at arbitrarily large Ri. 

In this paper we examine the behavior of stratified tur- 
bulence at large R values by means of numerical simu- 
lations of the Boussinesq Equations with forcing at the 
larger scales of the computational domain. Some sim- 
ulations of forced stratified turbulence have been per- 
formed previously by Herring and Metais 0, but the 
available resolution did not allow clear conclusions for 
flows at high R. Most experiments and simulations for 
stratified turbulence have been conducted on decaying 
turbulence, with many focused specifically on the transi- 
tion from isotropic to stratified turbulence after an initial 
high-energy excitation event pa, 0, El ■ Because of the 
large energy dissipation rate in both isotropic and strat- 
ified turbulence, this evolutionary path starts with large 
R and Fr and only briefly resides in a fully developed 
regime en route to small R and Fr; this situation there- 
fore provides a limited view of the equilibrium geophysi- 
cal regime. One of our main purposes here is to analyze 
the flow regimes of stratified turbulence in terms of R 
and Fr varied independently in a controlled fashion. 



II. DESCRIPTION OF THE CALCULATION 

A. Governing Equations 

A pseudo-spectral code is used to integrate the Boussi- 
nesq Equations on a triply periodic domain p^ : viz., 



d t u + u 
d t 6 + u • 
V • u 



Vu = -Vp+ 
V0 = ~N 2 u 
= 0. 



0e 3 
•e 3 + 



kA0, 



(2) 



In these equations, u is the three-dimensional velocity, p 
is the pressure divided by p„, k is the conductivity, and F 
is the imposed forcing. In Fourier space the equations for 
the Fourier components of velocity and the temperature, 



u(k) and O(k), are 

f (d t + vk 2 ) ui(k, t) - P l3 (k) 0(k, t) 
= -ifcjP in (k) J k+p+q=0 u„(p,i)u / (q,t)rf 3 p + F i (k,t), 

(d t + Kk 2 )e(k,t) + N 2 u 3 (k,t) 
„ = -i k„ J k+p+q=0 «n(p, i)0(q, t) d 3 p. 

(3) 

Pij = Sij — kikj/k 2 is the projection operator onto the 
plane orthogonal to u(k); 5y is the Kronecker tensor, 
repeated indices indicate summation; and i 2 — — 1 when 
i does not appear as an index. 



B. Energy Decomposition 

As a simple means of separating turbulence (rortices) 
and waves, as well as a horizontally uniform shear flow 
that is neither of these, we use Craya's decomposition 
for the incompressible velocity field in Fourier space into 
orthogonal components u„, u„,, and u s : 



u(k,t) = u w (k,t) + u U) (k,£) if k ft 7^ 0, 



u(k,t) = u«CM) 
where 



if k h = 0, 



u w (k,i) = 4> v (\s.,t) (fivQs), 
u w (k,t) =4> w Qs,t) <p w (k) 



(4) 



(5) 
(6) 



and 



^(k) = [(kxe 3 )]/|(kxe 3 )|, (7) 
<^(k) - [k x (k x e 3 ))]/|k x (k x e 3 )|. (8) 

k 3 = k • e 3 and k/ t = k — k 3 are the components of the 
wavenumber perpendicular and parallel to gravity. The 
components (4> V i4>w) were previously used by, e.g., Riley 
et al. and Lilly 01, and they are usually refered 
to "vortical" and "wave" components. The emergence of 
the "shear" component u s was emphasized by Smith and 
Waleffe j2£j. Associated with <f> v , (f> w and u s , we define 
the kinetic energy spectra by 



<M*0 = ^E&(p)&(p)' 

p 

M*) = \ J2 C(p)<Mp), 

P,Ph#0 



(9) 
(10) 

(11) 



p,Ph=0 



where the sum J] is done over a shell 
k - 1/2 < |p| < k + 1/2. We further define the "avail- 
able potential energy" spectrum by 



i ^ e*(p)6( P ) 



2 ^ iV 2 
p 



(12) 
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and the "total kinetic energy" spectrum by 

*(*) = 55>>) fi (p)- (13) 

p 

Total component energies (E v , E w , E s , E p , and E) are 
obtained by summing over all shells k. In addition to this 
decomposition, we define the "vertical kinetic energy" E z 
as half the area-averaged square of the vertical velocity, 
u • e3, and <S> Z as the vertical kinetic energy spectrum. 

C. Posing the Problem 
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Our analysis is based on a single numerical solution 
(designed after calculating many more exploratory so- 
lutions than we wish to admit). An anisotropic grid 
resolution is used to take advantage of the anisotropy 
of the flow that arises in response to the stable strat- 
ification (i.e., a small aspect ratio). The calculation 
is performed over a vertical fraction of a cubic domain 
by imposing a 2n/M vertical periodicity of the velocity 
(M = 8 in the present case). For a given number of grid 
points, this increases the achievable R value without loss 
of generality as long as the typical vertical scale is much 
smaller than the horizontal periodicity length. The level 
of stratification is controlled by adjusting the spatially 
uniform value of N in time. The forcing F is defined by 
F(k, t) dt = f3(k, t) u(k, t), with (3 chosen so that the dif- 
ference &f(k,t°) of the energy spectra before and after 
the forcing (i.e., the energy injection rate) is constant in 
time: 

<f/(M°) = i^{|u( P ,t) + ^F(p,<)| 2 } - $(M) 

p 

= (1+/3(M)) 2 *(M). (14) 
The coefficient /3(k,t) is obtained from eq. (|14|l as 



0(Js,t) 



l<f>(k,t) + $ f (k,t°) 



$(M) 



(15) 



We choose to force only the first vertical and horizon- 
tal modes (i.e., $ f (k,t°) ^ for k = k° = M and 
k = k° t = 1). In order to reach a high enough R value, a 
hyper-diffusion ((— l) p v p /^P h , with p=4 and a small coef- 
ficient of u p = 10~ 12 ) is added to the Newtonian dif- 
fusion in the horizontal direction. Several additional 
simulations at higher horizontal resolution demonstrate 
that this modification does not qualitatively affect the 
results presented here. Since most of the dissipation oc- 
curs in association with shear in the vertical direction 
(about 99.5% in the stable pancake regime; see below), 
we use ordinary Newtonian diffusion in this direction 
to have a clean definition of Taylor's Reynolds number, 
Rx = U rms \/(e/D) = 2VDE/e, where U rms = V2E 
is the RMS velocity, D is the enstrophy (i.e., vorticity 



FIG. 1: Experimental path for the primary simulation: the 
time evolution of Taylor's Reynolds number R\ (thick line), 
the resulting Brunt Vai'ssala frequency N (thin line), and the 
viscosity v (dashed line). N and v are adjusted to maintain 
the Froude number at the constant value of 0.08 and make 
R\ follow the indicated history. 



variance), A = ^/2E/D is the Taylor scale, and e the dis- 
sipation rate). The Prandtl number v/n is set to 1. The 
bulk Froude number is defined by Fr v — U rms /(HN), 
where H — 37r/4_E(^ fci k^ ^ v (ka)) is the typical verti- 
cal scale. 

The primary simulation is adjusted in time to follow 
a given experimental path in terms of R\ and Fr v by 
adjustment in time of N and v. The evolution of these 
two parameters is shown in Fig. ^ Fr v is maintained 
at a small value of 0.08, and R\ is changed from 200 
to 1000 by steps of 100 after integration periods of 100 
turnover times (i.e., T = L/U rms with L the integral 
scale defined by L = 3n/4E(J2 k fc _1 $(fc))). This period 
is approximately long enough to adjust to an equilibrium 
state for each R\. This experimental path is designed to 
expose a sequence of regime transitions with increasing 
R. 

The equations are solved on a 2ir x 2n x 2-k/M domain 
with a spatial grid size of 256 x 256 x 64 for R\ < 500 and 
256 x 256 x 128 for higher R\. The simulation appear to 
be slightly under-resolved for the highest Reynolds num- 
ber in the sense that the dissipation range is not extensive 
in <&(fc). The ratios of the highest resolved (fc™ = 96 and 
fc™ = 341) to the corresponding Kolmogorov wavenum- 
bers (k v h = (^/e^)- 1 / 4 and k% = (^/e^)- 1 / 4 ) are 1.1 
and 0.92, respectively, for R\ = 500, and they are 0.4 
and 0.2 for R\ = 1000. The under-resolution in the hor- 
izontal direction is compensated by the hyper-viscosity. 
In the vertical direction, the kinetic energy at the reso- 
lution scale is much smaller than at the largest scales in 
the dissipation range (by a factor of 10 -5 at the high- 
est R x ). The Ozmidoz scale (Lo = (e/N 3 ) 1 / 2 ~ 4.10~ 4 
for R\ — 700) is more than one order of magnitude 
smaller than the smallest resolved vertical scales. This 
ratio means that all the resolved scales are significantly 
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FIG. 2: Energy spectra with respect to horizontal and ver- 
tical wavenumbers averaged over more than 100 times at 
Rx = 200 for the primary simulation. 



FIG. 3: History of the energy components in the primary 
simulation. 



influenced by stratification. 

The simulations were made on one processor of a NEC 
SX-5. The code uses optimized ASL libraries of Fast 
Fourier Transforms and requires approximately 4 seconds 
per time step at the highest resolution. The time step is 
chosen to sufficiently resolve the fastest wave oscillations 
with period 2n/N. The simulation is integrated over a 
total of 360,000 time steps and takes about 200 hours. 



III. SOLUTION ANALYSIS 
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The experimental path for the primary solution is 
demonstrated in Fig. El as a function of the non- 
dimensional time, r = t/T, normalized by the turnover 
time T. The corresponding histories of the energy com- 
ponents are shown in Fig. [3J As R\ increases, we see a 
sequence of regime transitions. 

At moderate R\, less than ~ 400, the pancake motions 
are stable (i.e., the local Ri is everywhere large), and the 
vortical energy dominates all other energy components at 
all wavenumbers (Fig. [2J. The first transition occurs for 
R\ w 400, and it is evident in the significant growth of 
energy in the shear component, E s (Fig. |3J). The next 
transition, for R\ w 500, is evident in the intermittent 
occurrence of regions with small local values of Ri be- 
low the Kelvin- Hclmholtz critical inviscid stability value 
of 0.25 (Fig. 0). A further transition, for i? A « 700, 
is evident in local violations of the inviscid gravitational 
stability critical value of Ri = 0.0 (Fig. |3}. Finally, we 
see yet another transition for R\ w 900, evident in the 
growth of vertical kinetic energy, E z (Fig. |3J). Inter- 
estingly, throughout all these transitions, the principal 
measures of the internal- wave energy, E w and E p , show 
little change relative to the vortical-mode energy, E v . 
Since E v itself remains reasonably constant with time 
and its horizontal spectrum <fr v (kh) maintains a similar 
shape and magnitude at low wavenumbers, we conclude 



FIG. 4: Time evolution of R\ (solid line) and the volume 
fraction of the domain with local Ri < 0.25 (filled grey area) 
and with local Ri < (filled black area). There was no oc- 
currence of Ri < 0.25 for < r < 300. 



that pancake motions are indeed persistent throughout 
the R\ range we have been able to explore, even though 
the structure and intensity of the flow changes substan- 
tially at high wavenumbers, in the horizontally uniform 
vertical shear, and in the vertical velocity. 



A. Growth of the Vertical Shear Component 

In Fig. |3| the vortical energy is nearly steady over the 
entire simulation. The wave energy is more variable, but 
on average it is steady as R\ increases. However, the 
shear energy is a growing function of time. It repre- 
sents an inverse horizontal cascade of kinetic energy into 
kh = 0. The intensity of the inverse cascade of shear en- 
ergy is probably a function of the location of the energy 
peak in the horizontal direction; in the present case, the 
forcing is imposed at kh = 1, and a substantial part of 
the energy is transformed into pure vertical shear. To as- 
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FIG. 5: History of R\(t) (solid lines) and E 3 (t) for three sim- 
ulations. The primary simulation (with a dashed line for its 
Es) has an increasing Reynolds number by steps of 100 from 
200 to 1000, and the second and third simulations have con- 
stant Reynolds numbers of R\ = 300 (dash-dotted line) and 
R\ = 500 (dotted line), starting from the primary simulation 
at r = 200 and r = 400, respectively. 



FIG. 6: History of the global minimum value of local Ri 
for the primary simulation. The ordinate scale is split into 
three intervals to see both large and small values on the same 
plot. Whenever the minimum Ri drops below about 0.02, it 
continues to develop into overturning in the density profile 
(i.e., min[Ri] < 0); this first occurs around r ~ 400. 



sess the degree of equilibration for this inverse cascade, 
two additional simulations are made. Both start from 
the primary simulation and thereafter hold R\ constant 
for several hundred turnover times, but their starting 
(r, R\) values differ. Fig. [5] shows that E s does in- 
deed equilibrate over a period of less than 100 turnover 
times at a level that increases systematically with R\. 

This growth of shear kinetic energy has been seen 
previously when the Froude number is below a critical 
value [20| . For the alternative Froude number defined 
by Fr s — (e /c° 2 ) 1 / 3 /iV), our simulations have a value of 
approximately 0.025, more than an order of magnitude 
below the identified critical value of 0.42. In this previous 
study, the shear kinetic energy did not equilibrate even 
after more than 1000 turnover times. This may be due to 
its reliance on hyper-diffusion in all directions that exerts 
only a weak damping on the shear component. 



B. Onset of Overturning Motions 

Overturning occurs when an unstable shear layer rolls 
up, pulling high-density fluid above low-density fluid. 
This instability can occur in a fairly localized way. Its 
occurrence is detected with the local Ri defined in eq. 
|JT)|. and an indication of the overturning instability is a 
negative Ri value. Because of the Miles-Howard condi- 
tion, we can expect that this Kelvin-Helmholtz instabil- 
ity only occurs in regions where Ri is initially less than 
0.25. The history of the fraction of the domain with the 
local Ri below 0.25 is shown in Fig. 0] More events with 
Ri < 0.25 happen as R\ increases, although they evi- 
dently remain intermittent. The regime of overturning 
events first appears at R\ = 700 (r ~ 535). Even at the 



highest R value, only a small fraction of the domain is 
actively overturning at any time, less than 0.7%. The 
Probability Density Function (PDF) for Ri shows that 
most of the domain remains far away from overturning, 
with only a small tail in the PDF that extends to small 
and negative Ri values. Because of the intermittency, 
long averaging periods are required for stable statistics. 
Nevertheless, we can say that there is a well determined 
value of R for the onset of overturning events. To demon- 
strate this, the regime just before the first overturning 
event in the primary simulation (R\ = 500) is integrated 
over a longer period (250 turnover times), and no over- 
turning occurs despite several events with Ri < 0.25. To 
characterize more precisely the critical Ri for instability, 
we follow the history of the global minimum in local Ri 
(Fig. [fJJ). The minimum Ri value that does not immedi- 
ately lead to local overturning is Ri ~ 0.02. This value 
is a bit smaller but of the same order of magnitude than 
the values computed by Gage [|| for simple shear flows. 

The spatial distribution of small Ri events is organized 
into thin sheets with large vertical shear in the horizontal 
velocity. An example of a region with a negative vertical 
density gradient is shown Fig. [7| The vertical size of 
this region is very thin (a few grid cells) even though its 
horizontal size at this time (~ 27r/3) is not much smaller 
than to the domain size. An intense vertical velocity is 
associated with domain of negative vertical density gra- 
dient (Fig. [SJ, but the present simulation has only a 
marginally adequate resolution to expose the convective 
overturning events. 

Fig. shows the spectra for the energy components 
at t = 521, before any overturning occurs. The spectra 
with respect to horizontal wavenumber are very steep, 
with a slope close to k^ 5 for $„, and the shape does 
not vary much with R\ in this stably stratified regime 



FIG. 7: Instantaneous isolines of vertical density gradient at 
r = 741 (R\ = 900) for two perpendicular slices of the whole 
domain at z — cte (upper graph) and y — cte (lower graph) . 




FIG. 8: Vertical slice (same than in Fig. |7J of vertical ve- 
locity at t — 741 (R\ — 900) (upper graph) and a zoom of 
the region of intense vertical velocity with the projection of 
velocity vectors (lower graph) . 



without overturning. The onset of the overturning in- 
stability does not seem to affect significantly the overall 
level of wave and vortical energy (Fig. |3J). However, the 
overturning events are easily identifiable as a peak in the 
horizontal spectra at t = 741 (Fig. |3J). The peak is lo- 
cated at the horizontal scale Lh — 27r/50 that matches 
the typical vertical scale, L v . At this particular time, the 
instability is localized in a single region of the domain. 

The energy at the largest horizontal scales is not af- 
fected by the overturning instabilities that are localized 
in both space and time. In Fig. I1UI time-averaged 
vortical-energy spectra are compared for four values of 
R\. The horizontal spectra are very similar up to the 
typical scales of the overturnings. However, at both finer 
horizontal scales and at all vertical scales finer than the 
forcing scale, the spectrum amplitude increases system- 




FIG. 9: Energy spectra with respect to horizontal and vertical 
wavenumbers at r = 521 (R\ = 700) (upper graph) and r = 
741 (R\ — 900) (lower graph) for the primary simulation. 
The local energy peaks at r = 741 for 40 < < 60 are an 
indication of one or several overturning events(Min[ffi] > 0). 



atically with R\. At the constant, small Fr value in this 
simulation, the vertical spectrum slope becomes quite 
shallow as R\ increases. 



C. Growth of Large-Scale Vertical Motions 

The energy histories in Fig. [3] expose another tran- 
sition at an even larger R\ ~ 900, viz., the systematic 
growth of vertical kinetic energy E z . Inspection of the 
vertical energy spectrum reveals that the growth of E z 
after r = 700 occurs principally at k v = (Fig. Illfl . 
This mode of instability is reminiscent of the "negative- 
viscosity instability" observed in a Kolmogorov flow 0, 
further investigated by Dubrulle & Frisch Q with a 
multi-scale analysis. This analysis shows that a parallel 
flow with a small transverse scale develops a negative- 
viscosity instability to large-scale perturbations in the 
transverse direction when the viscosity becomes less than 
the RMS value of the streamfunction of the primary 
flow. In our simulation, it is difficult to test precisely 
this criterion of instability because the streamfunction is 
ill-defined. Crudely, we can expect an instability of this 
type when v < ^ s {d v ) 1 / 2 i v , where i v is the scale of the 
vertical shear and $ s {£v) is the shear energy at this scale. 
In our simulation this relation is satisfied on average for 
r > 700 if 2ir/£ v > 80. This scale is comparable in mag- 
nitude with the typical vertical wavenumber 2ir/L v ~ 50. 
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FIG. 10: Comparison of vortical (upper graph) and wave 
(lower graph) energy spectra with respect to horizontal and 
vertical wavenumbers for four different Reynolds number. 
Each spectrum is an average over more than 200 times for 
each R\. The horizontal spectra are very similar at large 
horizontal scales but differ at small scales (kh > 10). The 
situation is different in the vertical spectra, where the typical 
scale decreases with Rx- 
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FIG. 11: Evolution in time of the total vertical energy (thin 
line) and the vertical energy at the zero vertical mode (thick 
line) for the primary simulation (parameters of Fig. . The 
component of the zero vertical mode is several order of mag- 
nitude smaller than the total vertical energy for r < 700 and 
become the dominant mode for r > 800 



However, due to the complexity of the forced stratified 
flow, it is difficult to prove the nature of this instability 
pending more apt stability analyses. 



IV. SUMMARY AND DISCUSSION 

In our simulations of forced, equilibrium, stratified tur- 
bulence, we sec behaviors somewhat different from many 
previous studies of decaying stratified turbulence that 
were not able to sustain a large value of the Reynolds 
number, R. Most often the criteria for the occurrence of 
pancake vortices and suppression of overturning motions 
(i.e., Kelvin-Helmholtz and gravitational stability) have 
been linked to the stratification N but rarely to R. In- 
deed, we find that the stability of a solution is mainly 
controlled by two parameters with opposite effects on 
stability: increasing N (decreasing Fr) leads to a more 
stable solution and decreasing v (increasing R) has the 
opposite effect. For a fixed low value of Fr, we follow 
an experimental path of increasing R far enough to de- 
tect several regime transitions beyond the familiar one of 
stable pancake vortices. One transition is to the inter- 
mittent occurrence of regions with small or negative Ri. 
This refutes previous arguments [lii. ITfil ] that stratified 
turbulence remains stable with uniformly small local val- 
ues of Fr at large R and with uniformly cyclostrophic, 
hydrostatic diagnostic momentum balances. This tran- 
sition may plausibly be associated with inviscid Kelvin- 
Helmholtz and ensuing gravitational instabilities of the 
pancake vortices, although in our simulations the viscous 
effects on the unstable scales are significant. Neverthe- 
less, the pancake vortices continue to be the energeti- 
cally dominant component of the turbulence even up to 
the highest R values examined here, and visualizations of 
the large-scale potential vorticity field (not shown) show 
little change in spatial structure with R\. 

Two other transitions to large-scale motions other than 
pancake vortices do occur: a first one to growth of the 
shear kinetic energy at zero horizontal wavenumber and 
a second one to growth of vertical kinetic energy with 
its spectrum peak at zero vertical wavenumber for large 
R. The former has been seen previously in stratified 
turbulence |2jj, and the latter may be associated with 
negative- viscosity instability seen previously in unstrat- 
ified shear flow 0] ■ Each of these large-scale transitions 
may be interpreted as an inverse energy cascade. How- 
ever, they behavior is strongly constrained by the domain 
size in our simulations where the forcing is imposed at the 
gravest finite wavenumbers. We will explore beyond this 
limitation in future reports. 

In this paper we choose to focus on simulations at 
very small Froude number, and we are able to reach a 
Reynolds number high enough to destabilize the pancake 
vortices in several ways. This leads us to advance the fol- 
lowing proposition about the nature of equilibrium strat- 
ified turbulence: for any Froude number, no matter how 
small, there are Reynolds numbers large enough so that 
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a sequence of transitions to non-pancake motions will al- 
ways occur, and, conversely, for any Reynolds number, 
no matter how large, there are Froude numbers small 
enough so that these transitions are suppressed. Obvi- 
ously, this hypothesis warrants further testing, as do our 
provisional interpretations of the dynamical nature of the 
transitions. 



en Informatique Scientifique (IDRIS) and the IBM 
RS6000/SP of the Centre de Ressource Informatique of 
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